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A lattice model is presented for investigating the fluctuations in static granular materials under 
gravitationally induced stress. The model is similar in spirit to the scalar q-model of Coppersmith 
et al, but ensures balance of all components of forces and torques at each site. The geometric ran- 
domness in real granular materials is modeled by choosing random variables at each site, consistent 
with the assumption of cohesionless grains. Configurations of the model can be generated rapidly, 
allowing the statistical study of relatively large systems. For a 2D system with rough walls, the 
model generates configurations consistent with continuum theories for the average stresses (unlike 
W^" , the g-model) without requiring the assumption of a constitutive relation. For a 2D system with 

■ periodic boundary conditions, the model generates single-grain force distributions similar to those 

| obtained from the g-model with a singular distribution of q's. 
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The microscopic stress field in a static granular material has an extraordinarily complex structure. Viewed from 
the perspective of standard elasticity theory, the geometric disorder in the packing of grains gives rise to extremely 
OS j complicated boundary conditions on the stress equilibrium equations. In general, this disorder is rather difficult 
to characterize statistically and may even exhibit nontrivial correlations induced by the dynamical history of the 
material. 

At present, the overwhelming majority of attempts to model granular materials are formulated at the level of a 
continuous stress field, which is intended to represent a smoothed version of actual stresses, averaged over regions 
large compared to the grain size. In order to close the system of stress equilibrium equations, a constitutive relation 
■ must be assumed, such as the Mohr-Coulomb condition that the material is on the verge of yielding at everywhere 
within a plastic zone. While this condition jlj and others like it [^],[| have met with some success in many different 
situations, they rest ultimately on ad hoc assumptions about the connection between the microscopic and macroscopic 
stresses. 

One factor that could in principle pose fundamental difficulties for continuum theories for the average stress is that 

S fluctuations in the microscopic stresses may be quite large, perhaps large enough to invalidate typical assumptions 
about the scales over which the material can be modeled as a continuum. In several recent experiments that directly 
'■^J ■ image the stress field, stress chains (filamentary configurations of strongly stressed grains) have been observed with 
correlation lengths that are apparently limited only by the system size (though the systems have not been much larger 
than 30 grain diameters in the relevant dimension). For these reasons it is important to obtain some theoretical 
' understanding of what determines the size and spatial structure of the fluctuations. 
^ . Coppersmith et al. recently stimulated interest in a simplified statistical approach to stresses in granular materials 
with the introduction of the "q-model" , a lattice model whose configurations are intended to represent the way in 
which vertical forces are supported in a non-cohesive material. In such a model it is possible for very large fluctuations 
to arise, even on scales as small as a lattice constant. The constitutive assumptions of continuum models are replaced 
here by an ansatz concerning the form of the microscopic effects of geometric disorder in the material. The central 
result of Ref. || is that for infinitely wide layers (or materials subject to periodic boundary conditions in the horizontal 
directions) , fluctuations in the vertical forces supported by grains at a given depth are of the same magnitude as the 
average force supported at that depth. An elegant calculation shows that the probability distribution for the vertical 
force supported by a grain deep in the pile has an exponential tail, rather than Gaussian. 

For all its merits, the g-model has three serious flaws. First, it takes no account of the constraints imposed by 
horizontal force balance and torque balance on each grain, and thus does not contain the proper conservation laws at 
the microscopic level. This may constitute a flaw that leads to incorrect predictions, though it is also possible that the 
constraints in question do not affect the large scale behavior. Second, when studied in the silo geometry, the o-model 
yields predictions for the average stresses that dramatically disagree with classical theories and experiments. While 
quantitative agreement may not be expected given the crude representation of the boundary conditions that must be 
used in constructing gr-modcl configurations, the qualitative discrepancy is striking, as explained below. Finally, there 
is no clear procedure for connecting the lattice constant in the g-model to a physical length scale. 

In this paper, a new model is presented that explicitly incorporates the the relevant force and torque balance 
constraints into the lattice approach of the g-model and provides a natural connection between the lattice constant 
and the grain size. Stress distributions are then computed for the silo geometry with force-bearing walls and with 
periodic boundary conditions. It is shown that the new model gives much better agreement with previous theories 
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and experiments in the silo geometry, and thus appears to be a more reliable basis for investigating the subtleties of 
the stress fluctuations. Numerical results are then presented for the single grain weight distribution from the new 
model at large depths. The single-grain weight distributions are similar in form to those obtained in the g-model and 
therefore may be thought of as providing a firmer foundation for the g-model predictions. 

It is useful to make a conceptual distinction between the definition of the basic lattice with appropriate variables 
defined on it and the assumptions about those variables that are relevant for the study of non-cohesive granular 
materials. The geometric structure of the model is a square lattice with variables representing net normal forces, 
couples, and tangential forces on each edge. This structure will be referred to as the "net-lattice". Every possible 
stress field in any type of static medium can be mapped to a configuration on the net-lattice. (The mapping is 
many-to-one.) In order to model a non-cohesive granular material, several restrictions must be made on the values of 
the normal forces, couples, and tangential forces at each edge and an ansatz must be made for form of the disorder 
in the system. A particular model that incorporates these restrictions will be called the "a-model" , for reasons that 
will become apparent below. 

The paper is organized as follows: In Section Q the a-model is defined for 2D systems and the connections of the 
model parameters to physical parameters and constraints are discussed. In Section || results are presented for the case 
of narrow silos with rough walls in 2D and contrasted with results obtained from the g-model. In Section ^ results 
are presented for the case of wide 2D systems with periodic boundary conditions. Section |y| includes a discussion of 
some general issues and the generalization of the a-model to three dimensions. 



Consider an arbitrary (possibly highly inhomogeneous) material that is known to be in static stress equilibrium. 
For present purposes, the material is taken to be two-dimensional and the following lattice representation of the stress 
field is defined. (See Figure El.) A square lattice is constructed, with each cell representing a portion of the material. 
The length of each cell edge is 2r. A horizontal row of cells sharing vertices, such as those marked with a dot in 
Figure will be called a "layer" of cells. Cell ij is defined as the i th cell from the left in the j th layer from the top, 
with the left edge being defined as shown in Figure for even or odd j. Note that j increases downward, which is 
also defined as the positive y direction. 

Associated with each edge in the lattice are a normal force, a tangential force, and a couple. Taken together, these 
three quantities fully determine the force and torque exerted on one cell about its center by the other cell sharing that 
edge. The couple is necessary to represent the torque that can be applied to a cell even in the absence of tangential 
forces, due to the manner in which the normal force is distributed along the edge. The couple is here defined so that 
it does not include the contribution of the tangential force to the torque. The contribution from tangential forces is 
just equal to r times the net amplitude of the tangential force, regardless of how that force is distributed along the 
edge. 

Denote the amplitudes of the normal forces on the four edges of cell ij by N^', N^ 1 '^, n±\ n^'^'i the couples by 



Cj V, C% r, Ci r, Cg r; and the tangential force amplitudes by T^* , , t±\ 4 j as shown in Figure 0b. 



We will drop the superscript whenever this leads to no ambiguity. Both N and n refer to the amplitude of the net 
compressive force between the cells sharing an edge - negative values would indicate a net tensile force between them. 
The sign convention for tangential forces is chosen such that positive T or t always indicates a positive downward 
component of the force exerted by the cell with a higher center on the lower one. The sign convention for the couples 
is that the directions indicated in Figure 0b correspond to positive values of each of the Ci and CVs. This choice allows 
the right-left symmetry of the model to be immediately evident in the equations below. Note that the uppercase 
variables associated with a given cell are identically equal to the lowercase variables for the cell sharing the relevant 
edge. For example, n[ 1 '^ — n± (or n^~ lj 1 ') for j even (or odd). 

In terms of the canonically defined ^ microscopic stress tensor air), these forces and couples can be written as 
integrals over their respective edges. For example, letting ei and e2 be unit vectors in the y — x and y + x directions 
as shown in Figure 0b, and letting s run from — r to r along the edge: 



I. THE q-MODEL IN TWO DIMENSIONS 



A. The net-lattice 
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(2) 



2 




FIG. 1. Definition of the a-model. (a) The lattice of square cells. Shaded cells represent walls for the silo geometry and 
identified cells for periodic boundary conditions. Dashed edges are assumed to transmit no force in the silo geometry. The 
cells marked with a large dot constitute a single layer. The thick edges are used to compute the force on a layer to compare 
with the Janssen solution, (b) The variables used to describe the stress in a single cell and the unit vectors defining directions 
mentioned in the text. 
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where the integrals are over surfaces (lines in 2D) that may cut through grains and/or contacts between grains. 
Static equilibrium at the microscopic level requires [|| 

dk<Jik - P9% = 0, 



(4) 



where p is the local density and g is the gravitational acceleration. Requiring that the total force on a cell vanish, 
one finds: 



a l kdSk = mgt 



(5) 



where m = f pdV is the mass of the material in the cell and dSk is the fcth component of the outward normal to 
surface element S. Taking g to be in the positive y direction (downward) and performing the integral over the entire 
cell, the x and y components of Eq. (|| yield the following equations for vertical and horizontal force balance: 



nt + h + n 2 + t 2 = Nt + T x + N 2 + T 2 + V2mg, 
n% - h - n 2 + t 2 = Ni - Ti - N 2 + T 2 . 



Similarly, the vanishing of the total torque about the center of the cell requires 

<j) dS y x (a xy - (Jyy) + <j) dS x y {(jy X - o xx ) - ^ pgxdV = 0, 



(6) 
(7) 



which leads to 



h ~ ci - h + c 2 = -Tx -d+T 2 + C 2 + 2u, 



0) 
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where 2u — (y2/r) / dV pgx. (A torque about the center of the cell may be exerted by gravity because the center of 
mass of the material in the cell need not coincide with the geometric center.) 

The lattice and variables just described, together with the fundamental physical force and torque balance constraints, 
constitute the net-lattice. A "configuration" of the lattice denotes a set of values of , n 2 '^ c£' , c^'^ , t±\ 

4 <J) , and for all ij that satisfies Eqs. (§), (@), and (§) for every ij. Each configuration is a discrete 

representation of a possible stress field, and any possible stress field induces a configuration. Thus the net-lattice 
structure may be a useful tool for investigating stresses in a wide range of inhomogeneous materials. 



B. Definition of the a- model 



We now turn to the modeling of non-cohesive granular materials on the net-lattice. In the spirit of the g-model it is 
assumed that a typical stress distribution can be obtained by propagating the forces from top to bottom on the lattice. 
This assumption is made for the sole purpose of allowing the rapid construction of plausible force configurations. In 
the a-model, configurations are generated by descending through the lattice, solving Eqs. (JT^)- ( 13) for each cell. The 



process begins by specifying values of m and u for each cell and N\, N 2 , T±, T 2 , C±, and C 2 for the cells in the top 
layer. 

Given the values of the uppercase variables at a particular cell, the six lowercase variables must be determined. 
Since these variables are constrained by Eqs. (||), (m), and (^), the space of possible solutions is then three dimensional 
and can be parameterized by three real numbers ao, a\, and a 2 . The manner in which these three numbers are 
determined at each cell must reflect the physics of the material being modeled. 

It is convenient to make the definitions c\ = ai m, C2 = ot% ri2, and ao = T2 + C2 — ti + c\ + u. The force and 
torque balance equations for a single cell can then be written as 

m+h + n 2 +t 2 = N 1 +T 1 +N 2 + T 2 + V2mg, (10) 

n 1 -t 1 -n 2 + t 2 = N 1 -T 1 -N2 + T 2 , (11) 

ti - ain x + a = T 2 + C 2 + u, (12) 

t 2 - a 2 n 2 + a = T x + C\ - u, (13) 

where the torque equation has been split into two so that it may be expressed in a symmetric form. If m and u are 
given for every cell, then the microscopic true stress field induces a unique configuration of triples (ceo, a±, a 2 ) on the 
lattice. 

In a non-cohesive granular material consisting of convex grains, the grains can "push" on each other, but never 
"pull"; i.e., there can be no tensile forces in the material on scales larger than the grain size. This feature is 
incorporated into the a-model by imposing two restrictions: 

1. All normal force amplitudes n\ and n 2 must be positive; 

2. All ai's and a 2 s must lie in the interval (— 1, 1). 

Both restrictions are consequences of the assumption that there are no tensile forces anywhere along the cell edge. The 
maximal couple that can be produced by a given normal force amplitude corresponds to the case in which the entire 
normal force is applied at one corner of the cell, in which case a\ (or a 2 ) is ±1. The assumption is not rigorously 
valid since there can be tensile forces in the interior of a grain that is cut by the edge of a cell. The central hypothesis 
of the a-model, however, is that the geometric randomness in the grain positions can be replaced by the randomness 
in the choice of the solution for each cell. In this context, restrictions 1 and 2 are valid as long as the cell size is larger 
than or equal to the grain size. 

It is duly noted that the method of propagating forces down from the top of the system does not faithfully represent 
the physics in the following sense. The equations of stress equilibrium are elliptic equations whose solutions depend 
on the boundary conditions on all of the system's boundaries. For the propagation method employed in the a-model, 
the only cells that can be affected by a change in the distribution of forces on a given cell are the ones that lie in a 
downward opening cone with 45° edges, and in this respect the model is more appropriate to a hyperbolic system. 
The goal here, however, is to construct the simplest model consistent with the relevant stress balance constraints 
and capable of displaying nontrivial fluctuations. In addition, it should be noted that the classical approach to the 
computation of stresses via Mohr-Coulomb constitutive relations also transforms the problem into a hyperbolic one, 
and the computational method of propagating stresses downward from the top is routinely exploited in this context 
as well. 0] The development of rigorous elliptic methods for generating consistent configurations is an interesting and 
potentially important problem as well, but is beyond the scope of this work. 

The adoption of a hyperbolic method necessarily leads to a third restriction on the force amplitudes: 
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3. All tangential force amplitudes must be positive. 

In the absence of this restriction it is impossible to ensure that restriction 1 above can be met for every cell. That is, 
if negative tangential amplitudes are permitted to occur, then it can happen that for some cell deeper in the system 
all solutions to Eqs. (||), (0), and (||) have a negative value of either m or ni- On the other hand, it can be shown, 
using the method described in the following subsection, that restriction 3 guarantees a consistent solution for every 
cell. 

The final assumption defining the a-model is that geometric randomness in the structure of a granular material 
can be effectively modeled by assigning equal probability to all possible solutions to the force and torque balance 
equations at each cell. This requires defining a particular measure on the space of solutions. In the absence of any 
good reason to choose otherwise, it will be assumed that all values of the triple (ao, 0.1,0-2) that are consistent with 
restrictions 1-3 above are equally probable. Throughout the rest of this paper it will also be assumed, for simplicity, 
that all cells have the same m and that u = for every cell. 

This concludes the definition of the 2D a-model for non-cohesive granular materials consisting of convex grains. 
The next subsection presents a geometric picture of the solution space of the force and torque balance equations for 
a single cell. The picture provides useful insight into the need for restriction 3 and also reveals that the cell size in 
the model must be interpreted as the grain size. 

The g-model may be thought of as a rather drastic simplification of the a-model in which Eqs. ((7J) and (^) are 
ignored. The horizontal components of the forces can then be neglected and it becomes more convenient to work with 
variable Wi = (rij + tj)/ ' sqrtl and Wi = (Ni + Tj)/ 'sqrtl. Eq. (Q) can then be written as 

Wi = q(Wi + W 2 + mg) (14) 
w 2 = (l~q)(Wi+W 2 +mg). (15) 

In the g-model, q at each site is an independent random number chosen from some distribution with support only 
on the interval [0, 1]. The lack of correlation between g's at different sites is an important feature in the asymptotic 
analysis of Coppersmith et al. J6| Note that the a-model differs from the g-model in that the random variables at 
each cell cannot be chosen in advance, but only after the forces have been propagated down to the cell. If all of 
the Oi are fixed in advance, negative tangential and normal forces are quickly generated and the amplitudes diverge 
rapidly with increasing depth. Thus the "maximally random" a-model does contain correlations that are direct 
consequences of the additional force and torque balance constraints. In this sense it is related to attempts to derive 
force distributions by maximizing a statistical entropy associated with the possible ways of constructing microscopic 
configurations consistent with an imposed average stress field |fTo|| , but differs in that the average stress field in the 
a-model is not given in advance. 



C. Solving the force and torque balance equations in a single cell 

For the purposes of this subsection we will take g = 0. The effects of gravity can be included in a straightforward 
manner once the graphical method is understood. 

All of the information contained in the values of N, C, and T on one edge can be encoded in a single vector placed 
at some position along the edge. The vector itself represents the sum of the normal and tangential forces and the 
position is chosen such that the torque that would be produced about the center of the cell by the normal component 
of the force is equal to C. Thus the net effect of any configuration of stresses in the cell can be expressed graphically 
by drawing the four vectors, one on each edge. 

Now given the vectors on the upper edges of a cell we would like to determine all the possible ways of assigning 
vectors to the lower edges. One solution is immediately obvious and will be called the "direct solution" . Simply 
construct rays originating from the vectors on each edge, then take the lower vectors to be equal and opposite to 
the upper vectors and positioned where their respective rays intersect a lower edge. (See Figure ^j.) Each ray will 
be called a "ray of force" , though it is understood that the stresses represented are not really concentrated on the 
ray. The only complication that can arise is that both rays of force intersect the same lower edge. In this case, the 
vector assigned to this edge is the sum of the two, positioned such that the couple associated with it is the sum of 
the couples associated with the two. (See Figure g.) Note that the direct solution is guaranteed to exist if and only 
if each ray of force is guaranteed to intersect a lower edge, rather than exiting the cell through the other upper edge. 
This will be true whenever the normal and tangential amplitudes on the upper edge are both positive and is directly 
related to restriction 3 above. 

To construct solutions different from the direct solution, note that force and torque balance can be preserved by 
having a ray of force split into two rays at some arbitrary point in the cell, as shown in Figure pk. The force associated 
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FIG. 2. The direct solution for a single cell, (a) A straightforward case in which forces are propagated directly through the 
cell, (b) A case in which the propagated force rays intersect the same edge of the cell. The two vector forces must then be 
summed as shown to obtain the solution. 
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FIG. 3. splitting of force rays and tensile forces within a cell, (a) A vertex involving only positive force amplitudes, (b) A 
vertex involving a negative force amplitude, (c) A solution in a cell involving a negative force amplitude. The smaller dashed 
cells can be used to argue that real tensile forces must be associated with this configuration. 
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with each ray is directed along the ray and the force amplitudes are then fixed by requiring that the three forces 
sum to zero. Torque balance is guaranteed because none of the three forces generates any torque about the splitting 
point. Similarly, two rays of force that intersect can be merged into a single ray of force emanating from the point of 
intersection. 

Whenever all three rays of force intersecting at a single point lie in the same half-plane, the force amplitude 
associated with one of them must be negative. This situation is depicted in Figure [|b. The thin line in the figure 
indicates a negative amplitude. Any solution of Eqs. can be depicted as a network of force rays within the 

cell. An example involving a negative amplitude is shown in Figure ^p. 

It is possible for negative amplitudes to be an artifact of the summing of two forces that occur at different positions 
on the same edge, as would arise for the configuration shown in Figure ^p. There exist configurations, however, for 
which the negative amplitude can only be interpreted as corresponding to a real tensile force. An example is shown in 
Figure ||c. The forces applied to the upper edges of the cell result in no force being transmitted across edge E. The 
force transmitted across edge F balances both the forces and net torque produces by the forces above. There is no 
way to produce the configuration shown in a material that does not support tensile forces. To see this, consider the 
description of the configuration using cells of smaller scale (the dashed cells in Figure ||c. It can be shown that it is 
impossible to construct a network of force rays on the smaller cells that would be represented by the network shown 
on the large cell, without having tensile segments that cross the boundaries of some of the small cells. (The situation 
corresponds to that of a horizontal beam clamped at one end and supporting a load on the other, in which case the 
top portion of the beam is under tension.) 

Thus if all possible solutions of the force and torque balance equations occur with equal probability, the material 
within a cell must be assumed to support tensile forces and hence must be on the order of or smaller than the grain 
size. Together with the reasoning leading to restriction 1, this leads to the conclusion that the cell size in the a-model 
must be the grain size. 

In practice, the region in a-space which yields solutions that satisfy restrictions 1-3 is not easily identified. Restric- 
tion 2 immediately implies that a\ and a 2 must lie in the range (—1,1). Using Eqs. JlC|)-([L3"|) and the restrictions 
1-3, it can also be shown that cto must lie in the range 

- Ni - N 2 + Ci + C 2 - V2mg < 2a a < N x + N 2 + C x + C 2 + 2T X + 2T 2 - V2mg. (16) 

Thus the entire region of consistent solutions must be contained within a rectangular region in ce-space. To give 
equal weight to every solution consistent with the restrictions, a point is selected at random according to a uniform 
probability density throughout the rectangular box. The values of n\, n 2 , t%, and t 2 are then computed from Eqs. (|h})- 
(|l3|). If any of these quantities is negative, the solution is discarded and the process is repeated until a consistent 
solution is found. Since every point within the rectangular box has an equal probability of being chosen on every 
attempt, it follows that every point that corresponds to a consistent solution has an equal probability of being selected. 

In some cases, it may turn out that the set of consistent solutions occupies a very small region of the rectangular 
box so that the probability of finding a solution by random guessing is prohibitively small. The results described in 
this paper were obtained by imposing a cutoff K on the number of unsuccessful guesses. If the cutoff was reached for 
a particular cell, the direct solution for that cell was used. 

A final issue that must be addressed is how the gravitational force is to be distributed in the direct solution. Some 
care must be taken here to avoid generating more and more cells for which the direct solution must be used. The 
choice made for the simulations described in this paper is that the contributions due to the weight of the cell are 

m = V2mgii/(l + fi)(l + ai), (17) 
n 2 = V2mg/(l + (j,)(l + a 2 ), (18) 



ci = aim, (19) 

c 2 = a 2 n 2 , (20) 

h=c u (21) 

t 2 = c 2 , (22) 



where a\ and a 2 are each random numbers between 0.1 and 0.5, and \x = (N\ + T 2 + mg/v / 2)/(A 2 + T\ + mg/^/2). 
With this choice, the portion of the weight transmitted to the right (or left) increases when the cell is being pushed 
to the right (or left) from above. This is meant only to be plausible. It enters the computation of the stresses only 
when the direct solution is used, which can be made as infrequent as desired by raising the value of the cutoff K. 



7 



D. Boundary conditions 

The two sections below discuss the configurations generated under different sets of boundary conditions. In sec- 
tion ||, the focus is on the average stresses in a deep, narrow silo geometry, where walls can (and do) support stresses. 
The model then assumes the geometry shown in Figure [j]a with the shaded cells assumed to absorb all forces applied 
from above. The force and couple amplitudes on the dashed edges are taken to vanish identically. While there is no 
attempt to model the elasticity of the walls in a realistic way, significant differences between the a-model and the 
q-model will arise, leading to some useful insights. 

In section H, the goal is to explore the asymptotic behavior of the stresses in an infinitely wide and deep system. 
For this purpose periodic boundary conditions are employed in the horizontal direction. The two shaded cells in a 
layer in Figure 0a are identified and all cells are treated equivalently. The width of the system is always taken large 
enough so that the maximal force observed on any cell is small compared to the total force on its layer. 



II. AVERAGE STRESSES AND FLUCTUATIONS IN A NARROW SILO 

One classical approach to the computation of stresses in a granular material to assume that the material satisfies 
the Mohr-Coulomb criterion everywhere in space. Denoting the principal stresses as o\ and 02, the Mohr-Coulomb 
criterion in 2D reads 

■ =sm<£, (23) 

<J\ + <Ji 

where <fi is the internal friction angle characteristic of the material. A typical assumption for the boundary condition 
at the walls is a xy — taxi((f> w )a xx , where <p w is a parameter characterizing the friction between the material and the 
wall. 

Under these assumptions, the stress equilibrium equations can be solved numerically. Q The centerline stress a yy 
at depth y is well fit by the solution of Jansscn (which relies on additional simplifying assumptions) 



1 - e ay / R ) , (24) 



- _ Rpg 

a yy ~ „ 

with k = tan(0 w )(l + sin cj>)/(l — sin <j>). The Janssen solution assumes that a yy is independent of horizontal position, 
which does not turn out to be true in the full solution. ]5]|| Nevertheless, the deviations, which are of the order of 
15% from wall to centerline, do not alter the fundamental result that a yy /R is a function of y/R only. 

The behavior of the (/-model in this geometry is quite different. An ensemble average of configurations generated 
by the q-model yields a parabolic profile for the vertical forces supported. At large depths a steady state is attained 
in which the vertical stress near the wall must be proportional to the radius R of the silo, since the vertical force 
transmitted to the walls at each layer must equal the weight of the layer on average. The stress at the centerline, 
on the other hand, is proportional to R 2 . This behavior is accurately reproduced by an analytic calculation of the 
distribution of forces in the q-model when all <j's are assigned their mean value of 1/2. The analytic solution also 
indicates that the approach to the asymptotic profile occurs more slowly, with a decay constant of the order of 1 /R 2 
rather than l/R. 

Experimental results appear to confirm the scaling expectations of the classical analysis. The recent experiments 
of Clement et al., for example, show a linear scaling with system width in both the asymptotic average vertical stress 
and the characteristic depth of the approach to the asymptotic value. |ll]] 

Results for the average vertical force as a function of depth in the a-model are shown in Figure I Figure | shows 
data for three different silo widths, plotted in terms of the scaled variables suggested by the classical analysis. The 
average vertical force F v is here defined as the total vertical force transmitted across the corrugated surface marked 
with a thick line in Figure |l|a, divided by the number of cells in the surface, X: 

F v = ^y^^ ni + tl+n2+t2) (25) 

X 

Note that R — Xry/2. From Eqs. (0)-(||) the value of this quantity in the Janssen solution can be calculated. 
Accounting properly for the density of the material p — m/(2r) 2 and using the y = where the integer j indexes 

the layers, we find 

(26) 
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FIG. 4. Average vertical force on a layer in the silo geometry. The discrete symbols are obtained from simulation of the 
a-model. (See text for specification of parameters.) The thick dashed line is a fit to the Janssen solution, with one free 
parameter chosen to match the asymptotic value of the force at large depths. The thin line shows the behavior of the g-model 
with a uniform distribution of q's. 



where n is the combination of material and wall parameters defined above. 

Each data point in Figure |I] was obtained by averaging F v over 1000 configurations. In all cases, the cutoff K was 
taken to be 1000, resulting in the direct solution being used for approximately 3% of the cells. The data collapse 
obtained using the scaling suggested by the Janssen solution is quite good. The heavy dashed curve in Figure |l| is 
the classical prediction of Eq. (^6|) , where the single parameter a has been fit to the asymptotic value of F v at large 
depths. The same value of k = 2.15 was used for all data sets and the X was taken to be the average number of cells 
in two successive layers, not including the wall cells. It is perhaps worth noting that data from an experiment by 
Clement et al. show the same tendency to lie above the Janssen curve at small depths, [jflj The solid line in Figure |] 
is the prediction from the simplest version of the g-model, in which the distribution of g's is taken to be uniform over 
the unit interval, shown here for X = 39.5. The line shown was generated by simulation of the g-model, but agrees 
perfectly with the analytic solution. ]?J In that solution, F v scales at large depths like X 2 and there is no value of k 
that yields a satisfactory fit to the Janssen solution. 

It is clear that the inclusion of the proper force and torque balance constraints brings the a-model into much closer 
agreement with conventional expectations than the g-model. The reason for this appears to be that in the a-model 
strong stresses tend to propagate to the left or right with increasing depth, whereas in the g-model the vertical force 
simply diffuses. Figure || illustrates this difference with pictures of vertical force patterns obtained from both models 
for X = 39.5. 

The a-model clearly exhibits arching on the macroscopic level. The weight of the material is supported by the 
walls in the manner expected from continuum theories. It is also worth noting that the a-model exhibits well-formed 
arches on the "microscopic" level, an effect which will be illustrated more clearly below. In addition, the filamentary 
stress chains observed in the a-model are roughly reminiscent of experimental images produced under a variety of 
conditions (see, for example, Refs. [QH), and also with numerical solutions of stresses in random disk packings |l2| . 
Though quantitative comparison with these experiments and models is beyond the scope of this work, it does appear- 
that the model is capable of displaying plausible behavior both at the microscopic and macroscopic levels. 

Finally, it is instructive to consider the a-model in the absence of any randomness. If a\ and ai are set to 0.5 and 
the direct solution is employed for every cell, the model generates a smooth stress field in which the walls support 
almost no weight and the vertical force grows linearly with depth. This may be thought of as the correct result for 
a homogeneous solid with a vanishing Poisson ratio. The vertical stress generates no horizontal force on the walls, 
which therefore cannot bear any weight. Allowing a\ and 02 to vary between 0.1 and 0.5 but still employing the direct 
solution for every cell causes the walls to bear only a small fraction of the weight of each layer. 

This shows that the behavior illustrated in Figure ^ results directly from the randomness in the choice of solutions, 
not from the structure of the net-lattice and the "hyperbolic" method alone. No constitutive relation similar or even 
analogous to the Mohr-Coulomb condition of incipient yield has been put into the model, which makes the fact that 
the results agree reasonably well with the classical theory quite intriguing. 








FIG. 5. Typical configurations in the silo geometry. Darker cells indicate higher values of the vertical force supported. The 
bar indicates how the colors are specified on a linear scale. The width of the silo is X = 39.5. (Excluding the walls, there are 
39 cells in each odd layer and 40 in each even one. (a) The a-model. (b) The g-model. 
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III. FLUCTUATIONS IN A LARGE, PERIODIC BOX 



One of the primary motivations for developing the a-model is to see whether the predictions of the g-model survive 
the inclusion of physically realistic force and torque balance constraints. Three questions arc of particular importance: 

1. Does the probability distribution for stresses in an infinitely wide system approach a stationary limit at large 
depths? 

2. If so, does the limiting probability distribution have the same form as that of the g-model? 

3. Do the constraints induce any long range spatial correlations in the stresses? 

As a first step in answering these questions, distributions of the vertical forces were computed for a system of width 
500 at depths up to 450. The cutoff K was taken to be 1000, requiring direct solutions for approximately 3% of the 
cells. 

Figure || shows the P(w), the probability density for the vertical force supported by a single cell at various depths. 
The weight w is defined as the actual force supported, [n\ + ti + ri2 + ^2)/v / 2), divided by the depth of the layer. 
Part (a) shows P(w) for the deepest layer measured, along with the analytic result for the g-model with a uniform 
distribution of g's between and 1 which decays as exp(— 2w) at large w ||. 

In the a-model, it appears that P(w) decays as exp(— Aw) for large w with A roughly equal to 1.3, as shown in 
Figure ^|a. The dashed lines on the figure are guides to the eye, having slopes -2 and -1.3. Figure ^b indicates that 
there is a noticeable evolution of P(w) up to depths as large as 400, with the decay at large w becoming slower and 
the number of cells supporting almost no weight increasing with depth. 

The value of A obtained in the a-model can also be obtained from the q-model with a suitable choice for the 
distribution from which the g's are chosen. The open triangles in Figure |6|a were obtained from simulation of the 
g-model using a distribution of g's that included delta- function peaks at and 1; q was taken to be for 5.5% of the 
cells, 1 for 5.5%, and uniformly distributed between and 1 for the rest. The percentages were chosen in order to 
produce A = 1.3, as determined by the mean-field calculation described in the appendix. For comparison, Figure |?j 
shows the probability density for the values of q obtained in the a-model, where q is defined as the fraction of the 
vertical component of force on a cell that is transmitted to one of the cells in the next layer below, just as in the 
q- model. Note that the a-model does appear to generate singularities at q = and q = 1, but these are not delta 
functions. 

Figure |^ compares a measure of the spatial correlations in the a-model and the g-model. The correlation function 
(wiWi+j) c is plotted, where i indexes the horizontal position of a cell in a single layer. Though there is a discernible 
difference between the two models, it is clear that correlations decay rapidly, on the order of 10 cells. It may be argued 
that the q-model predictions should apply to the a-model at large w since the force and torque balance constraints 
do not appear to generate any long range correlations. Though there is as yet no analytic proof that the a-model 
weight distribution will indeed conform to the expected decay rate at very large w, the calculation in the appendix 
showing that this is the behavior expected in the g-model for an appropriately chosen q distribution, together with 
the fact that the q distribution obtained from the a-model is reasonably well approximated by a singular distribution 
of this type, strongly suggest the conclusion that the exponential decay with A w 1.3 will persist to arbitrarily large 
w. It must be noted, however, that the numerical data for the g-model appear to correspond to a slightly larger value 
of A. This may be due either to the influence of correlations not taken into account in the mean field caclulation, or 
possibly to the fact that the asymptotic decay rate emerges only for larger w or larger depths than were accessed in 
the simulations of Figure ^ Thus it is difficult to extract a more accurate value of A for the a-model from the data 
available at present. 

The maximally randomized a-model produces behavior more closely approximated by the q- model with 11% 0's 
and l's than by the maximally randomized q-model. Indeed, the two models yield rather similar spatial correlations 
as well as single site weight distributions. It is interesting to note that Radjai et al. reported an exponential decay 
in P(w) with A = 1.4 in numerical solutions for the stresses in 2D disk packings in squares of side length ~ 30 
disk diameters subjected to large external loading |1^| , not far from the value predicted by the a-model. In three 
dimensions, a similar singular q distribution was also found to agree best with dynamical simulations of spherical 
grains. || It is also worth emphasizing that in both the a-model and the q-model with an appropriate ^-distribution, 
there is significant evolution of P(w) for depths up to 450, even for small values of w. 

As shown above for the silo geometry, the a-model does provide a plausible picture of the macroscopic stress field. 
Unlike continuum theories, however, the a-model can also provide details on the scale of the grain size. Figure ^ 
shows a typical portion of a configuration at large depth for the a-model. Both the vertical and horizontal forces 
applied to each cell from above are shown. These images reveal that weight is supported by a network of arches 
with thickness on the order of the grain size. The appearance of such structures in a random model of this type is 
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FIG. 6. Probability distributions for vertical forces at various depths, (a) Comparison of the a-model with K = 1000 and the 
q-model with both a uniform distribution of g's and a distribution that is uniform except for sharp spikes at and 1 accounting 
for 11% of the total density. The data is from the layer at depth 450. The dashed lines are guides to the eye, with slopes -2 
and -1.3. (b) The a-model as in (a), but with data shown for several different depths. The inset shows the region near w = 0, 
where a clear evolution of P(w) with depth is observed. 
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FIG. 7. The q distribution generated by the a-model with K — 1000. Each data point represents the relative probability 
that the q of a given cell will fall in a bin of width 0.001 and the plot is normalized to correspond to rj(q) as defined in the 
appendix. The large circles mark the values for the bins centered on 0.0005 and 0.9995. The inset shows an expanded view of 
the singularity near q = 0. There is no measurable delta-function contribution at q = or q = 1. 
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FIG. 8. Correlations between weights on cells in one layer. The averaging was done over all cells at depth 200 in 9000 
configurations. Data is shown for the a-model with K = 1000 and K — 3 and the g-model with a q distribution that is uniform 
except for delta functions at and 1 accounting for 11% of the total density. 
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FIG. 9. Microscopic arching in the a-model. (a) Vertical forces supported by individual cells in the Q-model with K = 1000. 
Darker cells support larger vertical forces. The picture shown is a detail of a larger configuration, corresponding to a section 
the lower-most 70 layers in a system 240 layers deep. A clear central arch can be seen, together with several smaller arches, 
(b) Horizontal forces in the same region as (a). Darker cells are being pushed to the left by cells above them and lighter ones 
to the right. The arch apparent in (a) is seen here to have the expected structure of horizontal forces. 
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a nontrivial observation, as different ways of choosing (oto, cti, 02) can lead to substantially thicker chains and even 
nearly uniform distributions. 

Finally, a remark on the effect of changing K is in order. Changing K to 100 in the a-model generates direct 
solutions at 16% of the cells, but has little effect on the results described above. Changing K to order 1, however, 
creates a marked increase in the lengths and directions of the stress chains. A visual comparison of configurations 
obtained with K set to 1 or 1000 is shown in Figure |l(| 



IV. CONCLUSIONS 



The nrf-lattice is rich enough to describe the stress field in any material, with the scale of the cell size being 
completely arbitrary. For the case of non-cohesive granular materials with the cell size equated with the average grain 
size, however, simplifying assumptions can be made that lead to nontrivial predictions. 

The a-model studied in this paper includes particular choices of a few parameters that influence the details of the 
distributions. The latter category of choices has to do with the distribution of mass within each cell (the parameters u, 
mg), the precise form of the direct solution used when random attempts fail, and the assumption that all force incident 
on the wall cells is simply absorbed. Variations in how these choices are made might be expected to correspond to 
different choices for classical parameters such as the wall-material friction angle and the internal friction angle, which 
would be reflected in the value of k. 

Another parameter that can have a noticeable influence on the force configurations is the cutoff K that roughly 
determines how often the direct solution must be used. When K is small, the direct solution is used often and the 
details of how it is implemented can be important. For large K, however, the direct solution is used only in situations 
where the range of possible solutions is highly restricted anyway, so that all possible choices are actually quite close 
to the direct solution. For this reason, the configurations generated with K = 1000 are excellent approximations to 
the maximally randomized a-model. 

The present version of the model includes an assumption that permits configurations to be generated by propagating 
forces down from the top (restriction 3). This assumption is not rigorously justifiable, and may well be expected to 
fail in situations where strain effects are important. Savage has emphasized the importance of such effects, especially 
in the case of a free-standing pile. The boundary condition at the bottom of the pile, (the stiffness of the supporting 
substrate, for example) is known to be important in determining the stress field. fll3|| From the perspective taken in 
the current work, the question posed by the influence of the boundary conditions is how the boundary conditions 
affect the distribution of a's. Investigation of this issue might be possible if restriction 3 can be discarded and an 
algorithm developed for finding solutions consistent with appropriate boundary conditions on all sides of the lattice, 
including the bottom. In any case, the a-model is designed primarily to lend insight into microscopic and macroscopic 
fluctuations, not to investigate the details of how boundary conditions affect the average stress field. 

The solutions obtained from the a-model as constituted in this paper are sufficiently compatible with experiments 
on average stresses |llj and fluctuations |l4 ,[ll| to warrant further study. The a-model allows study of the qualitative 
features of the stresses at the grain size scale under the simplest physically consistent assumptions for the form of the 
geometric disorder. The effect of the disorder is taken to permit all possible solutions of the local force and torque 
balance equations with uniform probability in the solution space parameterized by ao, ai, and a2- Further work is 
needed to determine the sensitivity of the results to changs in the probability measure on this space. 

From the data presented in this paper, it appears that the 2D a-model predicts a weight distribution that decays as 
exp(— 1.3tu), consistent with g-model predictions if and only if an appropriate singular distribution of g's is used. For 
such a q distribution, the g-model also yields spatial correlations similar to the a-model. This may be taken both as 
an indication that the primary influence of the torque and horizontal force balance constraints is to select a particular 
form for the probability with which vertical force is transmitted between adjacent sites, and as a justification of the 
use of the g-model for understanding the basic structure of the stress fluctuations. 

In the silo geometry, the a-model achieves a possibly unexpected measure of success that is not obtainable by 
adjusting parameters in the q-model. The form of the average stresses generated by the a-model agrees rather well 
with experiments [jy]. This type of behavior arises also from the Mohr-Coulomb constitutive relation, which assumes 
that the material is at incipient yield everywhere. The a-model makes the rather different assumption that on the scale 
of the grain size the stress is as random as it can be without violating the fundamental conditions of stress equilibrium. 
The fact that this "works" suggests a conceptually new approach to the description of stress configurations in granular 
materials. 

Another intriguing connection of the a-model to recent work involves the explicit description of torques at the 
granular level. Experimental studies of the thickness of shear bands and also recent work on continuum models that 
include the dynamics of a field that characterizes the local rotation of the material, known as Cosserat continuum 
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FIG. 10. Comparison of the a-model with large and small K. The vertical force on each cell is shown is for periodic boundary 
conditions. At each layer the force is normalized by the layer depth, (a) K = 1000, resulting in use of the direct solution for 
3% of the cells, (b) K = 1, resulting in use of the direct solution for 80% of the cells. 
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models, have shown that shear bands tend to have characteristic widths on the order of 10 to 20 grain diameters. 
[p~6| -, P~7|| The occurrence of a similar length scale in a-model correlations suggests that the two approaches may be 
linked in a more profound way than has yet been understood. 

Generalization of the nci-lattice and a-model to three dimensions is straightforward but requires a substantially 
larger number of variables per cell. Using a face-centered cubic lattice oriented with the 111 axis on the vertical, one 
finds that there are 18 variables that must be computed for each cell. For each of the three downward-facing faces, 
one must find 

• a normal force, 

• two components of the tangential force, 

• two components of the couple (about the two axes that lie in the plane of the face) associated with the normal 
force, 

• and a third "torsional" couple (about the axis perpendicular to the face) determined by the spatial distribution 
of the tangential forces. 

The generalizations of Eqs. (^]), (^), and (^|) provide six constraints, one for each component of force and torque. The 
resulting 12-dimensional space of possible solutions can be parameterized by six a's relating the normal couples to the 
normal forces, three more analogous to the ag of the 2D model, and the three torsional couples. The high dimension of 
the solution space for a single cell makes the random guessing approach rather inefficient, and statistically significant 
data has not yet been obtained. 
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APPENDIX: MEAN FIELD CALCULATION OF A FOR THE g-MODEL 

Coppersmith et al. have derived several important results concerning the behavior of P(w) at large w in the g-model 
for various choices of the distribution of q's. [|| Define r/(q) as the probability that a given bond between cells will 
carry a fraction q of the vertical force exerted by the higher cell. Coppersmith et al. show that P(w) decays as 
exp(— Xw) with A = 2 in 2D for any 77(g) that has no singular contributions at q = (or q = 1). They also show, 
however, that different values of A can be obtained if such singularities are present. 

This appendix extends their mean field calculations to the case of distributions of the form 



V(Q) = ^(S(q)+ <%-!)) + (1 - 



(Al) 



in two dimensions, for which analytic results are obtainable. The term "mean field" here refers to the fact that 
all correlations between weights on adjacent sites are neglected. It is known that the mean field result is exact for 
certain special distributions, including the uniform one, and also that even in cases where it is not exactly correct, 
the deviations from it are small for large vertical forces. The calculation utilizes the Laplace transform formalism 
described in sections IID of Ref . || . Several results obtained there will be quoted here without explicit derivation. 

Let P(s) be the Laplace transform of the steady state weight distribution P(w) at large depths. P(s) satisfies 
equation 2.11 of Ref. M, reproduced here for the case of two dimensions only: 



P(s) 



dq 77(g) P(sq) 



(A2) 



In addition, normalization conditions imply that P(0) = 1 and P{s) — ► 1 — s as s — *■ 0. Substituting the desired form 
of rj yields 



P(s) 



-(1 + J P( S )) + (1 



1 -|2 

dq P(sq) 



" 



(A3) 



1G 



Defining R(s) = v P and changing variables in the integral, we have 



R(s) = £(i + j? 2 ^)) + (1-0)1 / dxR 2 (x). (A4) 



2 v w, • \ / s 

Multiplying by s and differentiating both sides gives 

dR _ l-R+(l-l)R 2 
ds s (1 - (W?) 

which can be solved for s in terms of i?, yielding 



(A5) 



bs= ~, r, (A6) 

(9 + (0- 2)i?)(2+e)/(2-e) ' 

where b is a constant of integration. Note that the normalization conditions on P imply R = 1 — s/2 in the vicinity 
of s = 0. Substituting this form for R in Eq. (|A6|), expanding about s = on the right hand side, and equating 
coefficients of the first order term yields 

b=^(26-2)-^^ 2 - e \ (A7) 

The inverse Laplace transform of P(s) will be proportional to exp(soiu) at large w, with So being the largest value 
of s for which P(s) has a singularity, which occurs wherever R(s) either has a pole or is of the form c + (s — so)" 
with the constant c ^ and the exponent v being non-integral. (If c — 0, then half-integral v also does not yield a 



singularity in R .) Although Eq. (A6) cannot be inverted to obtain R(s) explicitly, the position of the singularity in 



R{s) can be determined. First note that for < 9 < 1 Eq. (A6) implies that R(s) can diverge only at s = 0; since 
the coefficient of R in the denominator has magnitude greater than unity and the exponent is greater than one, the 
denominator must grow in magnitude faster than the numerator as \R\ goes to infinity. The possible divergence at 



s = arises only because we multiplied by s to obtain Eq. (A5). (As mentioned above, R is known to approach 1 at 



s = 0.) Next note that dR/ds must diverge at any value of s for which R = 1/9, and only at those points, as is evident 



from Eq. ( A5) and the fact that R itself does not diverge. Finally, repeated differentiation of Eq. (A5) reveals that 



higher derivatives of R cannot diverge at any p oint where the first derivative doe s not diverge. Thus the singularity 



in R can be located by setting R = 1/8 in Eq. (A6) and combining with Eq. (A7) to determine s. The result is that 
the singularity occurs at 

so approaches as 9 approaches 0, which may be expected given that 9 = corresponds to the critical distribution for 
which the q- model exhibits power law decay in P(w). Q Moreover, so approaches -2 as 9 approaches 1, indicating a 
smooth convergence to the result derived for the uniform distribution. (The case of 9 = must be treated separately, 
however, and it is seen that R develops a pole at -2.) Finally, the exponent v at the singularity may be obtained from 
lim^jyg \n(d 2 R/ds 2 / \n(dR/ds) = (y — 2)/{v — 1), which yields v = 1/2. Thus there is a single singularity in R(s) 
and in the vicinity of the singularity we have R = (1/9) + y/Js — sq), with sq given above. This result is consistent 
with the claim in Ref. || that P(s) has a square root singularity for any r)(q) having a delta- function component at 
q = 0. 

In order to compare to the n ume rical results for the a-model, it is useful to find the value of 9 that produces a 
decay with A = 1.3. From Eq. ( |A8| ) one sees that so ~ —1.30. . . is produced by 9 = 0.11, which is the reason that 
this value was chosen for plotting in Figures 6 and 8. 
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